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In the adiabatic post-Newtonian (PN) approximation, the phase evolution of gravitational waves (GWs) from 
inspiralling compact binaries in quasicircular orbits is computed by equating the change in binding energy 
with the GW flux. This energy balance equation can be solved in different ways, which result in multiple 
approximants of the PN waveforms. Due to the poor convergence of the PN expansion, these approximants tend 
to differ from each other during the late inspiral. Which of these approximants should be chosen as templates 
for detection and parameter estimation of GWs from inspiraling compact binaries is not obvious. In this paper, 
we present estimates of the effective higher order (beyond the currently available 4PN and 3.5PN) non-spinning 
terms in the PN expansion of the binding energy and the GW flux that minimize the difference of multiple 
PN approximants (TaylorTl, TaylorT2, TaylorT4, TaylorF2) with effective one body waveforms calibrated to 
numerical relativity (EOBNR). We show that PN approximants constructed using the effective higher order 
terms show significantly better agreement (as compared to 3.5PN) with the inspiral part of the EOBNR. Eor 
non-spinning binaries with component masses mi 2 6 [1.4Mq, ISMq], most of the approximants have a match 
(faithfulness) of better than 99% with both EOBNR and each other. 


I. INTRODUCTION AND SUMMARY 

One of the biggest scientific enterprises of recent times, the 
quest for the direct detection of gravitational waves (GWs), is 
expected to achieve its first success in the near future. Some of 
the second-generation of interferometric GW detectors [1,2] 
will start operating later this year and are expected to achieve 
their design sensitivity over the next few years [3]. Estimates 
of the astrophysical rates of candidate GW sources, in partic¬ 
ular coalescing compact binaries, predict that first detections 
are within the reach of these observatories [4]. 

GW signals from coalescing compact binaries, buried in the 
noisy detector data, are to be detected by cross-correlating the 
data with theoretical templates of expected signals. These 
theoretical GW templates are computed by solving the field 
equations of General Relativity for the two body systems. 
The two-body problem in General Relativity has no analyti¬ 
cal exact solution. Hence the construction of GW templates 
either requires approximation techniques that help to tackle 
the problem analytically, or large-scale numerical computa¬ 
tions for solving the problem exactly. In the early stages of 
the inspiral of the compact binaries, where the orbit can be ap¬ 
proximated as an adiabatic sequence of quasi-circular orbits, 
the GW templates can be computed using the post-Newtonian 
(PN) approximation to General Relativity [5]. However, the 
modeling of the dynamics of the system and the GW signals 
from the late inspiral and merger stages requires large-scale 
numerical relativity simulations [6]. (See, e.g., the discussion 
in [7] about the region in which the PN description of the in¬ 
spiral is valid.) 

For “low-mass” binaries (total mass < 12Mo), the GW sig¬ 
nal observed by ground-based interferometric detectors will 
consist almost entirely of the inspiral portion of the wave¬ 
form [8, 9], which could, in principle, be modeled accurately 
by the adiabatic PN approximation. In this approximation, 
the phase evolution of the binary’s quasi-circular orbit (and 
hence the gravitational waveform) is computed by equating 
the loss of the orbital binding energy with the energy flux of 
the GWs. This energy balance equation can be solved in dif¬ 
ferent ways, which result in multiple approximants of the PN 
waveforms. Due to the poor convergence of the PN expansion, 
these approximants tend to differ from each other during the 
late inspiral stage. Which of these approximants should be 
used as templates for detection and parameter estimation of 


GWs from inspiraling compact binaries is not obvious. Luck¬ 
ily, for the case of non-spinning binaries, most of the standard 
PN approximants are shown to be effectual [10] for the pur¬ 
pose of GW detection; however they are not faithful to the 
actual signals for accurate estimation of the source parame¬ 
ters [9]. For the case of highly spinning binaries, on the other 
hand, the currently available PN approximants fail to be even 
effectual for GW detection using advanced detectors [11]. 

Motivated by this issue in GW data analysis, we seek to 
estimate effective higher order terms in the PN expansion of 
the binding energy and GW flux such that multiple, if not all, 
PN approximants have close agreement with a fiducial “exact” 
waveform family. We do so by fitting the PN approximants of 
an appropriate dynamical quantity [which we choose to be the 
evolution of the PN time t{u) as a function of the PN expansion 
parameter v\ with that computed from a fiducial exact wave¬ 
form over a range of mass ratios. This fitting is done over the 
putative inspiral regime — frequencies less than that of the 
Schwarzschild innermost stable circular orbit (c = 1 / Vfi). As 
the fiducial exact waveform family, we choose the effective 
one body waveforms calibrated to numerical relativity (EOB¬ 
NR v2) [12]. We estimate two effective higher order terms — 
or “pseudo-PN” (pPN) terms — in the binding energy (log- 
independent and log-dependent terms at 5pPN order), and five 
effective higher order terms in the GW flux (log-independent 
term at 4pPN, and log-independent and log-dependent terms 
at 4.5pPN, and 5pPN order). These are shown in Fig. 1. Cur¬ 
rently we restrict ourselves to the case of non-spinning bina¬ 
ries. 

We show that multiple waveform approximants (TaylorTl, 
TaylorT2, TaylorF2) generated using the 5pPN accurate en¬ 
ergy and flux functions show excellent agreement (faithful¬ 
ness ^ 0.99 - 0.999) with the fiducial exact waveform fam¬ 
ily (EOBNRv2) over the whole “low-mass” parameter space 
mi ,2 e [1.4Mo, ISMq] (see Figs. 3 and 4). The TaylorT4 ap- 
proximant shows very good agreement (faithfulness ^ 0.99) 
over most, but not all, of the parameter space. For the majority 
of the cases, the faithfulness of the pPN approximants are sig¬ 
nificantly better than that at 3.5PN order (see Fig. 3). Since the 
other adiabatic PN approximants are all basically variants of 
TaylorTl, TaylorT2, and TaylorT4, we expect these results to 
hold for other approximants also. These pPN coefficients can 
be readily applied for searches for GWs from non-spinning 
low-mass binaries. Work is ongoing to extend this method to 
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the case of spinning binaries. 

We stress that the higher order terms that we estimate are 
effective terms, in the sense that they capture the effects of a 
large number of actual higher order terms in the PN expan¬ 
sion. Hence they will be different from the actual terms at a 
given PN order. This is in contrast with the work that has been 
done in the extreme mass-ratio limit, where one is able to de¬ 
termine true PN coefficients in the flux and binding energy (in 
addition to spin precession and tidal effects), since one can 
work to very high precision (up to thousands of digits) and 
at very large radii (up to 10^° times the Schwarzschild radius 
of the central black hole), allowing one to easily disentangle 
the individual PN coefficients, and to extract their analytical 
expressions (see [13, 14] for application of these methods to 
the binding energy). We also note that there is related work 
by Huerta et al. [15], who similarly fit for effective PN coef¬ 
ficients using BOB waveforms in the context of intermediate- 
mass-ratio inspirals. 

The rest of this paper is organized as follows: Section II 
briefly reviews the adiabatic energy-balance equations, as well 
as the computation of different PN approximants of the GW 
phase evolution. Sec. Ill introduces the pseudo-PN coeffi¬ 
cients, describes our method for determining them, and then 
evaluates their performance by computing mismatches be¬ 
tween EOBNR waveforms and PN waveforms. The conclud¬ 
ing section IV summarizes the results and suggests possible 
avenues for future work. We give various ancillary techni¬ 
cal results as appendices. Throughout this paper we use ge- 
ometrized units: G - c - 1. 


II. POST-NEWTONIAN WAVEFORMS IN THE 
ADIABATIC APPROXIMATION 


A compact binary revolving about its center of mass radi¬ 
ates away orbital binding energy via the emission of GWs. 
This radiation in turn causes the separation R between the 
component masses to shrink, and the orbit decays. For a large 
portion of the binary’s inspiral, the rate of shrinkage of R (or, 
equivalently, the rate of increase of angular speed) is negligi¬ 
ble over the duration of an orbit. During this adiabatic regime, 
the rate of loss of binding energy £ may be equated with the 
GW flux ff radiated.* For a quasicircular binary, this relation 
alone results in a set of coupled ordinary differential equations 
for the orbital evolution, called the phasing formula [17]: 

dt ~ &'{vy dt ~ M’ ^ ’ 


where M is the total mass of the binary. The orbital binding 
energy &{v) and the GW flux ffiv) can be computed as PN 
expansions in terms of a gauge-invariant velocity parameter 
V (see, e.g., [5] for a review). These equations, which describe 
the time evolution of the orbital phase (p and the velocity pa¬ 
rameter V, may be written in integral form as: 


f(y) = fref + 



, dv, 

r(v) 


(2.2a) 


^ Actually, one can only equate the GW flux to the system’s mechanical en¬ 
ergy loss up to the Schott term. However, this term vanishes for exactly 
circular orbits and is thus small for quasicircular ones when one is in the 
adiabatic regime; see, e.g., the calculation of the Schott contribution in the 
EOB framework and associated discussion in [16]. 


<f(v) = (frsf + — I v^——dv. (2.2b) 

M T{v) 

These integral expressions describe the time and phase evolu¬ 
tion of the orbit as a function of the expansion parameter v. 
Here, iijef is a reference value of v while fjef and represent 
the time and phase of the orbit at c = Cref- 

The phasing formula can be solved in a number of differ¬ 
ent methods which are perturbatively equivalent, in the sense 
that they all are accurate to a given PN order. They can be 
generally classified into three approaches: 

1. Compute the energy and flux functions in Eq. (2.1) up 
to a given PN order, evaluate the ratio 'F{v)fS'{v) nu¬ 
merically, and solve the ordinary differential equations 
using an appropriate numerical method. The resulting 
time-domain approximant is known as “TaylorTl” [17]. 
One can also compute an equivalent frequency domain 
approximant, “TaylorFl” making use of the stationary 
phase approximation [17]. 

2. Re-expand the ratio 'F(v)IS'{v) in Eq. (2.1) as a power 
series and truncate it at the respective PN order. The re¬ 
sulting time-domain approximant obtained by solving 
the ordinary differential equations numerically is called 
“TaylorT4” [18]. Frequency domain equivalents of Tay- 
lorT4, such as “TaylorF4” and “TaylorR2F4” can be 
computed via the stationary phase approximation [11]. 

3. Re-expand the ratio S'{u)lff{v) in the integral form of 
the phasing formula (2.2) as a power series and truncate 
it at the respective PN order. This allows us to evaluate 
the integrals analytically resulting in a parametric repre¬ 
sentation in terms of t(u) and (fiiv), from which one can 
obtain ip{t). This is known as the “TaylorT2” approx¬ 
imant [17]. Other time-domain approximants making 
use of this re-expansion include “TaylorT3” [17] and 
“TaylorTS” [19]. Based on this re-expansion, one can 
also compute a frequency domain approximant, “Tay- 
lorF2,” making use of the stationary phase approxima¬ 
tion [17]. 

The energy function &(v) and the flux function ffiv) are 
known only to a finite PN order: For the case of non-spinning 
binaries, S(v) is known to 4PN order [20, 21], and ffiv) to 
3.5PN order [22, 23]. (See also [5] for a review of these com¬ 
putations.) The approximants therefore increasingly diverge 
from each other as the inspiral progresses and v gradually in¬ 
creases to a considerable fraction of the speed of light. This 
divergence is not acute during early-inspiral (u <s: 1), but can 
become considerable during late inspiral [11, 18]. An obvious 
conundrum that arises is the choice of approximant to make 
to better model the “true” GW waveforms of Nature during 
late inspiral, leading to merger. Waveforms that more accu¬ 
rately represent GWs from inspiraling compact binaries ex¬ 
ist, such as those produced by large-scale numerical-relativity 
(NR) computations. In principle one could use these wave¬ 
forms and do away with PN waveforms. However, this is not 
feasible in practice, given the enormous computational cost of 
NR computations. EOBNR waveforms are an excellent sub¬ 
stitute for NR waveforms (at least in the non-spinning case). 
There is unfortunately still a non-trivial computational cost 
associated with producing EOBNR waveforms, in view of the 
fact that millions of them will be required to construct a tern- 
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plate bank suitable for GW observations.^ In the next section, 
we propose to introduce effective higher order terms in the PN 
expansion of the binding energy and GW flux as a method to 
overcome this. 


III. INTRODUCING PSEUDO-PN TERMS IN THE 
ENERGY AND ELUX 


Circumventing the onerous task of computing higher order 
PN terms analytically (which would still not have necessar¬ 
ily improved the accuracy of the approximants), we introduce 
a single set of higher order effective PN coefficients in the 
energy and flux of the binary to simultaneously improve the 
agreement (as quantified hy faithfulness [ 10 ]) between wave¬ 
forms produced by multiple approximants, and their EOBNR 
counterparts. 

The PN formalism expresses the binding energy and flux 
as expansions in powers of v. So far, binding energy PN 
coefficients up to eighth order in v beyond the Newtonian 
order (4PN) have been determined; furthermore, for non¬ 
spinning binaries, coefficients at half-integer PN orders are 
known to vanish through 4.5PN, but start to be nonzero at 
5.5PN [13, 28, 29]. On the other hand, PN flux coefficients 
have been computed up to 3.5PN, and unlike for the binding 
energy, there are nonzero coefficients at all half integer PN 
orders starting from 1.5PN. 

We propose determining the pPN coefficients for &(u) and 
ff (v) up to 5pPN order, starting from the lowest order at which 
PN coefficients have thus far not been determined. We assume 
the following ansatz for our pPN coefficients, guided by the 
form of known PN terms, some of which include quantities 
proportional to In v (known from the test-mass limit of the en¬ 
ergy flux—see, e.g., [30] —and the first-order self force results 
for the binding energy given in, e.g., [28]): 


£5pPn(c) = --Mr]v^ 


^£‘ kv '' + (£10 + eJ'o 1 ” 


k=0 


(3.1a) 


32 


,2..10 


'^SpPn)^) - V 
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k=0 k=% 


(3.1b) 


A. Determining the pseudo-PN coefficients 

As mentioned earlier, waveforms generated by large scale 
NR computations are the fiducial waveforms of choice, to 
which we could calibrate our PN waveforms and determine 
the pPN coefficients. However, NR waveforms spanning the 
long inspiral are computationally expensive, and are thus not 
yet available for binaries with arbitrary mass ratios. In order 
to lift this restriction, we use instead the EOBNRv2 waveform 
model [12] as implemented in the LALSimulation package, 
part of the LALSuite software library [31]. 

We broadly considered two ways by which we could es¬ 
timate our pPN terms. In the spirit of TaylorTl, the first 
method we considered attempts to extract pPN coefficients 
by fitting the ratio - ffiu) f&'(ii) (keeping the energy and flux 
to 5pPN order) to the velocity derivative dv/dt computed us¬ 
ing EOBNR, over a range of u spanning the late inspiral, say, 
V e [0.2, uisco]. where uisco = 6 “'^^ is the velocity at the in¬ 
nermost stable circular orbit (ISCO) for a test particle orbiting 
a Schwarzschild black hole. 

As expected, the agreement between TaylorTl and EOBNR 
improves noticeably. Eor example, the faithfulness for mass 
combination 1.4Mo - 1.4Mo increased from about 0.95 at 
3.5PN to about 0.999 at 5pPN. However, the pPN coefficients 
so determined were significantly larger (by three to four or¬ 
ders of magnitude) than the average size of PN coefficients.^ 
Other approximants computed using these pPN terms exhib¬ 
ited poor faithfulness with the EOBNR waveforms. Including 
the known terms from the test-particle limit in the energy and 
flux (to 22PN, using the exact energy and the flux from [30]) 
did not tame these pPN coefficients - they still remained un¬ 
desirably large. 

The second method, which, after some tuning, yielded 
promising results, computes t(u) from Eq. (2.2) using Tay- 
lorT2, up to 5pPN order. The coefficients of t{v) at orders 
beyond 3.5PN are functions of the pPN quantities (Eio, E\q, 
^8.9.10, £ 3 ^, 9 , jq]. (We chose this number of coefficients since 
it is the smallest number that gave us the level of agreement 
we desired in the final matches, and also includes all the ex¬ 
pected terms in both the energy and flux at a given PN or¬ 
der; we have not experimented extensively with adding fur¬ 
ther coefficients.) Instead of fitting the t{v) we obtain from 
these coefficients to its corresponding EOBNR analogue, and 
thus estimating the energy and flux pPN coefficients directly, 
we choose instead to define a functional form for the TaylorT2 
pPN terms as follows: 


Here, 77 := 11111712 !fvf' is the symmetric mass ratio of the bi¬ 
nary, where m\ and m 2 are the individual masses, ffk are 
the known PN coefficients, and £ 10 , £^q, F^, are the pPN 
and pPN-log terms to be determined by calibrating PN quan¬ 
tities such as t(v) or (f{v) given in Eqs. (2.2) to their EOBNR 
counterparts. Just like ordinary PN coefficients, we expect 
pPN terms to vary smoothly with the symmetric mass ratio 77 , 
ideally as low order polynomials in 77 . 


f5pPN(l') - 


5M 

25677778 




\k=0 



where 


08 = In (; + 6^ In^ u. 


(3.2) 


(3.3a) 


^ This issue is partly solved by the development of surrogate models of 
EOBNR waveforms [24] making use of reduced-order modeling tech¬ 
niques [25-27]. However, the construction of these reduced-order models 
requires the generation of tens of thousands of EOBNR waveforms, and 
then finding an orthonormal basis for them. This also incurs a significant 
computational and memory cost. 


^ One can gauge the order of magnitude of PN coefficients by looking at the 
test particle limit, where these coefficients are known to high orders. One 
finds that the «PN coefficient increases in size roughly as 3”'^^, due to the 
divergence at the light ring (see, e.g.. Fig. 3 in [14]), so that at the pPN 
orders we are considering, we expect flux coefficients on the order of 10^ 
to 10^ and binding energy coefficients on the order of 10^ to 10^; see the 
explicit expressions for the test particle flux in [30] and the general form of 
the test particle binding energy in, e.g., Eq. (3.3) of [10]. 
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FIG. 1: The estimated pseudo-PN coefficients (red points) plotted against the symmetric mass ratio rj of the binary, along with the polynomial 
fits to these data points given by Eqs. (3.6). (We do not plot Fg here, since we set it to the test-particle value for all rj.) 


0 ^ = 0 *^ H- 6 ^ In y, A:e{9,10). (3.3b) 

This ansatz closely parallels the actual TaylorT2 expansion of 
t(v) (see [9] for the expression through 3.5PN and [32] for the 
expression to higher orders in the extreme mass-ratio limit); 
as a result, to the orders we consider, the 6 coefficients are all 
linear in terms of the pPN coefficients of the energy and flux. 
The higher order terms in the t{v) expansion can be written 
in terms of the higher order PN coefficients in the energy and 
flux (to six significant digits) as: 

= 2730.98 - 1015.307/ - 225.8907/^ - 45.02017/^ 

- 44.01047/"^ -I- 8Fs, (3.4a) 

= -386.268 - 1564.247/ -h 4^^, (3.4b) 

09 = -12572.5 -H 15468.47/ -h 10627.87/^ -h 2328.76?/^ -h 8Fg 

- 8 ^, (3.4c) 

0^ = 3278.27 H- 8 ^-^, (3.4d) 

010 = 15242.5 - 5463.777/ -H 3469.497/^ -i- 863.8287/^ 

- 402.8247/^ - 33.94467/^ H- (23.6905 22.6667t])Fs 

- (11.8452-H 11.33337/)^^-24£io + 10£^o+47^10 

- 2F\o, (3.4e) 

05^0 = -1951.32 - 6083.577/ - 4606.307/^ -h (23.6905 

-I- 22.66617])F^ - 24E\o + 4/^^ . (3.4f) 

The next step involves fitting our ansatz for the fspPNCf) 
given in Eq. (3.2) to its EOBNRv2 counterpart fEOBCf) for a 
range of mass ratios. We compute v = from the dom¬ 

inant quadrupole {{ - m — 2) mode 1)22 of the EOBNRv2 
waveforms, where the orbital frequency w is computed in the 
following way: 

1 dwoo 

W = - — 7 —, (022 = arg(t)22)- ( 3 . 5 ) 

2 at 

Here we compute the derivative using second-order centered 
finite differencing of the numerical data. We then fit the result¬ 
ing data for fEOB(o) to the analytical expression fspEN)!*) using 
a least-squared minimization algorithm, to obtain numerical 


values of the pPN coefficients 0^, 0^ and 0jr^ for different mass 
ratios. Here we use a least-squares fit (i.e., minimize the 
norm of the difference of the functions) not only because it 
is a standard fitting procedure, but also because we are pri¬ 
marily interested in improving the matches between the func¬ 
tions, which are fairly closely related to the norm."^ Ad¬ 
ditionally, as discussed in Appendix A, the pPN coefficients 
obtained from a least square fit in t{u) automatically improve 
the least square residual of (fiv). In performing these fits, we 
have the freedom to set the EOB and pPN t{v)s to be equal to 
each other at a reference value of v by adding a constant: We 
choose to do this at y = 0 . 2 . 

While an obvious approach would be to fit for all these co¬ 
efficients over the entire range of y values we are considering, 
we found that we obtained better results (possibly closer to the 
actual PN coefficients) if we used a more involved method. In 
this method, we realize that the lowest-order coefficients we 
consider (i.e., Eg and Eg at 4pPN) will be dominant at small 
y, so it makes sense to only fit for those coefficients over a 
restricted range. However, we do not know this range a pri¬ 
ori, so we let the upper limit of the interval over which we 
fit (the transition velocity yg) vary and minimize the residual 
over the entire range of y we consider [Umin, Umax], to determine 
Ug (which we take to be at least slightly larger than Umin)- We 
then subtract off the contribution of the 4pPN terms and move 
on to the 4.5pPN terms, where we apply the same method, fit¬ 
ting over the interval [Umin, Ug], requiring that y'g < Ug < Umax- 
Eor the 5pPN terms, we fit over the entire interval, since we 
want to improve the overall agreement as much as possible, 
and are not adding on any more terms. This iterative method 
was inspired by the iterative method for obtaining the PN co¬ 
efficients in the linear-in-7/ portion of the binding energy from 


'* Indeed, by Plancherel’s theorem, in the case of white noise, maximizing the 
match is equivalent to minimizing the norm of the difference between 
the functions, since in this case the match is just the inner product be¬ 
tween the (L^) normalized functions. In fact, under appropriate simplifying 
assumptions, one can obtain an explicit lower bound on the match between 
the two waveforms in terms of the norm of the difference of their phases, 

given in Appendix A 2. 
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FIG. 2: The difference between t(v) as computed via EOBNR and TaylorT2 (both rescaled by the Newtonian value of f), at orders from 3.5PN 
to 5pPN, for three different mass ratios. The improvement in the agreement between EOBNR and TaylorT2 (due to the addition of pPN terms) 
becomes more dramatic with increasing mass ratio. 


a high-precision self-force calculation used in [14]. The fits 
we quote are done using Umin = 0.2 and Umax = Osco- 

The fit-and-minimize-over-transition-velocities method we 
just described yields some unwanted jumpiness in the pPN co¬ 
efficients, the largest components of which can be attributed to 
a similar jumpiness in the dependence of the transition veloc¬ 
ities on rj. In order to reduce this jumpiness, we thus fix the 
values of the transition velocities to Ug - v’^ - 0.355, which 
are very close to the values for these quantities returned by the 
above procedure for an equal-mass system (and a number of 
other mass ratios; indeed, almost all the u^s given by the min¬ 
imization are extremely close to 0.355). We find that these 
fixed transition velocities give final residuals at 5pPN that are 
quite similar to the residuals obtained from the fits that allow 
the transition velocities to vary, while reducing much of the 
unwanted jumpiness in the pPN coefficients. 

We also ended up not including a variable coefficient of 
In^(u) at 4pPN in the fit, since when we did so, we found that 
the values assigned to the correction to the test mass value in 
the 4pPN In(Li) term in the energy flux by the fit were small and 
quite jumpy. We thus fixed the coefficient of In^(y) at 4pPN by 
taking the 4pPN In(Li) term to just have its test-mass value for 
all j] and then went back and performed the fit-and-minimize- 
over-transition-velocities procedure again to arrive at the final 
transition velocity values quoted above. 

We display the residuals we obtain with this method as a 
function of v in Fig. 2, which plots the difference between t as 
computed via EOBNR and the (p)PN t{v) values (both scaled 
by the Newtonian value of f) for (p)PN orders from 3.5PN 
to 5pPN. The plot shows that the residuals get progressively 
smaller with increasing PN order beyond 3.5PN, a clear in¬ 
dication that the pPN coefficients systematically improve the 
agreement between TaylorT2 and EOBNR. This improvement 
in the residuals should thus translate into an improvement in 
matches, as discussed above. 

Having determined the pPN coefficients 
0 g, 0g^, 0g, 0g, 010, 6 * 10 , we can extract the energy and flux pPN 
coefficients via Eqs. (3.4). It is obvious, upon inspecting 
these equations, that the flux coefficients at 4 and 4.5pPN 
(viz., Eg,Eg,Eg,Eg) are uniquely determined. However, 
the coefficients Eio,Eio’ Eio,Ej'g are not fixed, and may in 
fact take any value provided they satisfy their constraint 
equations (3.4e) and (3.4f). Exploiting this freedom in the 
5pPN energy and flux parameters, we endeavor to set them in 
a way so as to improve the agreement between the TaylorTl 
and EOBNR waveforms. To that end, we fit the 5pPN flux- 


to-energy ratio T(u)IS'{v) to EOBNR’s dvidt, after fixing 
the pPN coefficients at 4 and 4.5pPN to those determined via 
the method described above. We use the constraint equations 
Eqs. (3.4e) and (3.4f) to write the 5pPN flux coefficients as a 
function of the 5pPN energy coefficients, and then perform a 
two-parameter minimization over the energy coefficients. 

We compute the pPN terms in the energy and flux func¬ 
tion making use of a set of 21 EOBNRv2 waveforms, with 
symmetric mass ratio values evenly spanning the interval 
rj e [0.05,0.25], as shown in Eig. 1. Eitting low order polyno¬ 
mials in T] to the numerically computed pPN coefficients (and 
including the test-particle value to which we set Eg for all 77 , 
for completeness), we write these pPN coefficients (to six sig¬ 
nificant digits) as: 

Eg = -229.100 - 934.582/7 - 861.481/7^ 

E^ = 52.7431, 

Eg = 1146.18 - 2743.15/7 - 22150.5/7^ -H 64309.3/7^ 

F\ = -421.553 - 595.925/7 - 15568.9/7^ h- 48222.8/7^ 

Eio = 11051.9 h-21160.8/7-215292/7^-1- 18123.0/7^ 

= 15725.6 - 4275.47/7 - 157235/7^ - 279859/7^ 

£10 = 1966.08 -H 6752.30/7 - 56757.4/7^ -h 34764.9/7^ 

E\o = 2565.84 h- 2946.63/7 - 44485.0/7^ - 21789.4/7^ 

(3.6) 

We plot these fits, along with the original pPN coefficients at 
discrete 77 in Eig. 1. We see that these coefficients are generally 
of a reasonable size (and the original 0 coefficients—which 
we do not show explicitly — are generally of similar size to 
the test-particle TaylorT2 coefficients at these orders, given 
in [32]), except for the binding energy coefficients, which are 
a bit large. 

Numerical methods employed: In the first part of the esti¬ 
mation of the pPN coefficients, where we fit the EOBNR t{v) 
using a TaylorT2-like ansatz, we minimize the residuals over 
the transition velocities using the Nelder-Mead downhill sim¬ 
plex algorithm implemented in Scipy’s [33] optimize.minimize 
function, and determine the pPN coefficients at each order 
using the Levenberg-Marquardt non-linear least square algo¬ 
rithm implemented in Scipy’s optimize.curve _fit function. In 
the second part, where we determine the energy and flux pPN 
coefficients at 5pPN by minimizing the difference of a Tay¬ 
lorTl expression for dv/dt with that computed from EOBNR, 
we compute the EOBNR dv/dt using second-order centered 
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finite differencing. In order to alleviate the modulations in 
our numerically computed EOBNR dvidt (likely coming from 
the residual eccentricity in the EOBNR waveform due to dif¬ 
ficulty in setting exactly quasi-circular initial conditions) we 
have generated all the EOBNR waveforms we consider start¬ 
ing from a velocity of y = 0.15, while we only fit for v > 0.2. 
Eurthermore, we smoothen the numerically computed v{t) by 
employing Scikit-Learn’s [34, 35] implementation of a ma¬ 
chine learning algorithm called isotonic regression. 


B. Testing the performance of pPN coefficients 

To evaluate the performance of our method of obtaining 
pPN coefficients in the energy and flux functions, we compute 
the mismatch (1 - faithfulness) between the different approx- 
imants computed using these pPN coefficients and the cor¬ 
responding EOBNR waveform. Eollowing the discussion in 
Sec. II, we consider three approximants — TaylorTl, Tay- 
lorT4, TaylorT2 — which correspond to the three different 
ways of treating the ratio of the energy and flux functions ap¬ 
pearing in the phasing formula. In addition, we also consider 
the frequency-domain TaylorE2 approximant. Although the 
behavior of the TaylorE2 approximant is expected to be very 
similar to that of TaylorT2, we explicitly consider the former 
because it is the approximant that is most widely used in GW 
data analysis.^ We indeed see that the behavior of this approx¬ 
imant is very similar to that of TaylorT2. 

We evaluate the mismatches over the frequency range 
[10 Hz, /isco], where 10 Hz corresponds to the lower cut-off 
of Advanced LIGO’s expected frequency band, and /isco is 
the dominant quadrupole mode GW frequency associated with 
the ISCO in the Schwarzschild geometry. The mismatches are 
computed assuming the “high-power, zero-detuning” noise 
power spectral density of Advanced LIGO [36]. We consider 
a set of component masses oti -m 2 consisting of standard fidu¬ 
cial systems 1 AMq -1 AMq, IOMq - 1 AMq, IOMq - IOMq, as 
well as three others, to better sample the space of total masses 
and mass ratios. (Note that the lower bound of the frequency 
band over which the mismatches are computed is smaller than 
the lower bound, v = 0.2, of the range of v over which the fits 
were conducted to compute the pPN coefficients for all the 
mass combinations we considered. The smallest u that is in 
band is y ^ 0.08 for the 1 AMq - 1 AMq system.) Eor each of 
these systems, we compute the mismatch as a function of the 
(pseudo) PN order, starting from OPN and working our way up 
to 5pPN.® The results are summarized in Pig. 3. We see that 
the pPN terms reduce the mismatch of all the approximants 
with EOBNR for almost all the cases considered here; in most 
cases the mismatches have been reduced to < 10“^. 

The only instance where the pPN terms have slightly wors¬ 
ened the mismatch is in the case of TaylorT4 approximant for 
comparable-mass binaries with large masses (q ^ I, M > 
IOMq). Given that none of the pPN coefficients were de¬ 
termined by improving the agreement between TaylorT4 and 


^ Note that TaylorF2 is based on the stationary phase approximation, which 
obtains additional corrections starting at 5PN, as discussed in [32], which 
we do not include here. 

® For the sake of consistency, we set the amplitude order equal to the PN 
order whenever possible, otherwise defaulting to the maximum amplitude 
order available in LALSimulation [31], which is 3PN for everything except 
TaylorF2, which only uses the Newtonian amplitude. 


EOBNR, it is not surprising that the pPN coefficients are rel¬ 
atively less effective at improving their agreement. [In fact, 
TaylorT4 and TaylorT2 are, in some sense, maximally dif¬ 
ferent from each other, since the former is based on the re¬ 
expansion of ;F(y)/£(y), while the latter is based on the re¬ 
expansion of its inverse.] However, it is worthwhile to no¬ 
tice that the pPN coefficients still improve the matches at high 
mass ratios. Eurthermore, for comparable mass systems, the 
mismatch at 5pPN, though marginally worse than at 3.5PN, 
still remains below 10“^. 

Having looked at a set of discrete mass combinations, we 
further test the performance of our pPN terms by evaluating 
the mismatches over a continuous two-dimensional region of 
the nil -m 2 parameter space, with mi 2 6 [2Mo, 15Mo]. Com¬ 
paring the mismatches at 3.5PN with those at 5pPN (Pig. 4), 
we find, for three of the four approximants considered, a 
consistent reduction in mismatches at 5pPN (as compared to 
mismatches at 3.5PN) over the entire mi - m 2 region dis¬ 
played. Eor TaylorT4, the improvement in the mismatches 
occurs mostly for higher mass-ratio systems, while it does not 
become worse than ~ 10~^ at comparable masses. 

Also of interest is the match between approximants. Given 
that the PN expansions of the binding energy £(y) and flux 
T (y) are known only to a limited PN order, the approximants 
are not expected to converge towards each other, a fact that 
becomes manifestly evident during late inspiral. Since our 
pPN coefficients improve the match of various approximants 
with EOBNR (often considerably), we thus expect that they 
will also help the approximants better converge towards each 
other. Pigure 5 shows the mismatch between different pairs of 
approximants as a function of the (pseudo) PN order. The fig¬ 
ure demonstrates a significant improvement in the agreement 
between TaylorTl and TaylorT2. Since the pPN coefficients 
do not reduce the mismatch between TaylorT4 and EOBNR 
quite as significantly as for the rest of the approximants, it is 
not surprising that the effect of the pPN coefficients in amelio¬ 
rating the difference between TaylorT4 and the rest is not as 
marked; nevertheless, we still find improvement in the agree¬ 
ment for a good fraction of the mass combinations considered. 
Note that since TaylorT2 and TaylorP2 are closely related ap¬ 
proximants, and given the similarity between the TaylorT2- 
EOBNR and TaylorP2-EOBNR mismatches, we did not find 
it necessary to include mismatches involving TaylorP2 in the 
plots. As expected, we found the latter to be qualitatively sim¬ 
ilar to mismatches involving TaylorT2. 


C. Understanding the resnits 

As discussed in Appendix B, it is not guaranteed that pPN 
coefficients that improve one approximant will necessarily im¬ 
prove another one, since different approximants (e.g., Tay- 
lorT2 and TaylorT4) are expanding quantities with different 
analytic structure. In particular, poles in one case become 
zeros in the other, so the expansions can have, in principle, 
very different radii of convergence. However, in Appendix B, 
we also obtained sufficient conditions for the pPN coefficients 
in the flux that improve TaylorT2 to also improve TaylorTl, 
which is easier to ensure than requiring them to also improve 
TaylorT4, since the pPN version of TaylorTl can be thought of 
as a rational (including logarithms) representation of the pPN 
version of TaylorT2 (similar to a Pade approximant). Here we 
saw that having small pPN terms at each order (as we obtain 
with the current method, but had difficulty obtaining with the 
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FIG. 3: Mismatch with EOBNR as a function PN order, for approximants TaylorTl, TaylorT2, TaylorT4, and TaylorF2. The legend shows the 
masses in Mq of the binaries considered. For most of the approximants and mass combinations considered, the mismatch at 5pPN is not only 
better than the mismatch at 3.5PN, but is also lower than 10“^. 



FIG. 4: Contour plots of the mismatch of different approximants with EOBNR over the mass range nii 2 e ISMq]. The top panels 

correspond to 3.5PN order and the bottom panel to 5pPN order. There is a consistent improvement in the mismatches at 5pPN, as compared to 
3.5PN, for TaylorTl, TaylorT2, and TaylorE2, over the entire low-mass region examined here. The mismatches for TaylorT4 improve mostly 
for high-mass ratios; nevertheless, for most the low-mass region, the mismatches are near or below 10“^. 


methods that we previously tried) is sufficient for tfiis to be 
the case. Indeed, we find that if we extract the pPN flux co¬ 
efficients tfirough 5pPN from the TaylorT2 9 coefficients, set¬ 
ting tfie 5pPN energy coefficients to zero, then these pPN flux 
coefficients already improve the agreement between TaylorTl 
and EOBNRv2. Since we then go on to fit for the 5pPN en¬ 
ergy coefficients in TaylorTl, it is not surprising that the final 
5pPN TaylorTl has excellent agreement with EOBNRv2. 


IV. CONCLUSIONS AND OUTLOOK 

Accurate models of GWs from inspiralling compact bina¬ 
ries are a crucial component in the searches for these signals. 
Given that the construction of template banks will require mil¬ 
lions of waveforms to adequately cover regions of the bina¬ 
ries’ parameter space, minimizing computational costs is an 
important challenge. Numerical relativity is able (in princi¬ 
ple) to supply accurate waveforms that give the predictions 
of general relativity. However, such waveforms are extremely 
expensive to generate, particularly considering the many cy¬ 


cles necessary to completely cover the band of a GW detec¬ 
tor [37]. EOBNR waveforms, which calibrate an analytical 
model to numerical relativity, reproduce the numerical rela¬ 
tivity results quite well over large portions of the parameter 
space and can be generated in a small fraction of the time 
it takes to produce a numerical relativity waveform. How¬ 
ever, they still require solving a nontrivial system of ordinary 
differential equations, and thus are still too expensive to use 
to construct template banks directly in most cases. On the 
other hand, PN waveforms are quite inexpensive to generate. 
However, they are unreliable during late inspiral where they 
begin to diverge from each other as well as from their more 
accurate NR or EOBNR counterparts. To ameliorate this is¬ 
sue in a manner that does not significantly increase compu¬ 
tational costs, we introduce effective higher order PN terms 
(which we call “pseudo-PN” terms) in the energy and flux 
functions of the binary designed to improve the agreement be¬ 
tween EOBNR (our fiducial waveform of choice) and multiple 
PN approximants. 

We use a procedure described in Sec. Ill, which, in essence, 
fits a PN-like ansatz for t(v) and dii/dt to the corresponding 
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FIG. 5: Mismatch between approximants as a function PN order. The legend shows the fiducial mass-combinations (in Mq) that we consider. 
There is a significant reduction in the mismatches at 5pPN (which are consistently below 10“^), as compared to 3.5PN, between TaylorTl and 
TaylorT2. The mismatches involving TaylorT4 also drop for all cases considered. 


quantities computed from EOBNR, to estimate the effective 
higher order terms in the energy and flux. We compute co¬ 
efficients for an evenly spaced set of symmetric mass ratios 
rj € [0.05,0.25], and find that the pPN coefficients seem to 
vary as low order polynomials in rj (Fig. 1). To evaluate the 
performance of the estimated ty-dependent pPN coefficients, 
we compare the mismatches of approximants TaylorTl, Tay- 
lorT2, TaylorT4, and TaylorF2 evaluated at 5pPN order with 
our fiducial “exact” waveform FOBNR (see Figs. 3 and 4). 
We find that for approximants TaylorTl, TaylorT2, and Tay- 
lorF2, the mismatches with FOBNR over mass-space mp 2 6 
[1.4Mo, 15Mo] are not only smaller than those at 3.5PN, but 
are around or below 10“^, often touching 10“^. The pPN coef¬ 
ficients do not significantly improve the performance of Tay- 
lorT4, but they do nonetheless reduce the mismatch for high 
mass ratios and maintain a mismatch of close to 10“^ for a 
significant portion of the parameter space considered. Mis¬ 
matches between approximants also indicate a marked im¬ 
provement at 5pPN, as compared to 3.5PN, for most of the 
cases considered. Furthermore, a large fraction of these yield 
mismatches of below 10“^ at 5pPN (Fig. 5). 

Based on the encouraging results produced by our pPN co¬ 
efficients in boosting the performance of the approximants we 
have considered, we anticipate similar enhancements in other 
approximants closely related to ours (see the discussion in 
Sec. II). We believe therefore that the pPN coefficients will 
afford some flexibility in the choice of approximants, were 
they to be used as part of the pipelines for the detection and 
parameter estimation of GWs from non-spinning sources. 

Of course, the most useful application of our method for 
GW data analysis would be to extend to the case of spin¬ 
ning binaries. We do not foresee this to be particularly chal¬ 
lenging in principle (at least in the case of binaries with non- 
precessing spins), and plan to pursue this as an extension of 
the work presented here. In addition, as successful as FOBNR 
is at reproducing NR waveforms, it may not always paral¬ 
lel full-scale NR in terms of accuracy and reliability. This 
is particularly the case when spinning systems are considered 
[38, 39]. Computing pPN coefficients with NR waveforms 
may therefore be a worthwhile endeavor, given that a number 
of long and accurate NR waveforms spanning the full param¬ 
eter space of interest are becoming available (see, e.g., [40]). 
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Appendix A: Bounds on the residual for <p(v) from the residual 

for t{v) 

Here we wish to show that, if we minimize the residual 
between the pPN and EOBNR versions of t{v), it automati¬ 
cally guarantees a small residual for ifiv). We will first mea¬ 
sure the residual using its maximum, i.e., using the norm 
ll/IU max„e[ 0 ,„„„] \f(v)\ for / : [0, y^ax] ^ where we can 
replace the usual supremum by a maximum since /’s domain 
is compact (and do not have to consider the essential supre¬ 
mum since we are only applying it to continuous functions). 
Later we will also use the norms for p 6 [1,2], defined 

for a general p e [l,oo[ by ||/||/, := \f(vWdvf^‘’ (we 

will use a lower bound of y = iimin in the discussion of the 
matches). The properties of these norms, including the stan¬ 
dard inequalities we use here, are discussed in most books on 
real and/or functional analysis (e.g., [41]). 
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1 . L” bounds 


We now assume that we have found fpPN(y) such that ||Af||oo 
is small (where A denotes the difference of the pPN and BOB 
quantities), and we want to bound ||A^||oo, where ip'{v) - 
(v^lM)t'{v) [from Eq. (2.1)]. We can thus integrate directly 
and write 


M||Ai/j||oo = max 


f 




(Al) 


where tc and ipc denote the time and phase at coalescence, re¬ 
spectively, and (/, g )2 f(t)g{t)dt is the standard in- 

ner product for real functions on the interval [fmin, Wx]- Now, 
we can write the waveforms as /ji(f) = A[v{ty) cos 
where v{t) is the same for both waveforms, by assumption, 
so \\hk \\2 — ||A o v\\ 2 l V 2 , where we have written cos^ifk - 
{cosltpic + l )/2 and neglected the contribution of the oscil¬ 
latory cos2(fiii piece. If we now write A(f ;= tpi - ip 2 , so 
cos (fi cos ip 2 - [cos Aif + cos((,oi -H ^ 2 ) 1 / 2 , and again neglect 
the oscillatory cos(i,tJi -H ^ 2 ) term, we have 


We now integrate by parts to remove the derivative on At and 
then apply the triangle inequality (and triangle inequality for 
integrals) to obtain 


M\\Aip\\^ < max 

l«[0,Omax. 




|Af(y)| + 3 I |Af(i;)| dv 


(A2) 


We are able to discard the boundary term at y = 0 in the in¬ 
tegration by parts: While t{v) diverges as as y \ 0, both 
the pPN and BOB f(y)s reduce to the PN t(v) for small v and 
the leading 4pPN term we are adding is 6 >[ln^(y)] for y \ 0, 
so we know that v^At(v) vanishes at least as fast as In^(y) as 
y \ 0. We now note that we can bound the maximum of the 
sum by the sum of the maxima of each term, and each of these 
maxima can be bounded from above by y^j^^HAfUoo separately. 
Bor the hrst term, this is trivial, and for the second, almost 
so, after noting that one can bound the integral from above by 
its value with |Af(y)| —> ||Af||oo, where the latter is a constant 
we can pull out of the integral (or, alternatively, applying the 
Holder inequality in the form Il/gl1 1 < II/II 1 II 0 IIOO). Our hnal 
result is thus (recalling that we take u^nx - yisco = 6 ^'^^) 

M||A(^|U < 2yLxllAf|U = ^IIAflU - 0.136||Af|U. (A3) 

This inequality is surely not optimal, though it certainly suf- 
hces for our purposes. Note that we can replace the lower 
bound of y = 0 by a lower bound of y = yo > 0 , as we use 
in our hts, with no change, since we perform the hts setting 
Af(yo) = 0. 


2 . LJ bounds 

We are also interested in bounds involving the norm, 
since these are related more directly to matches, as mentioned 
in Sec. Ill A. However, going from the At and Aip which we 
have been considering to bounds on anything involving the 
waveform is not entirely straightforward, since one has to in¬ 
vert f(y) to obtain the waveform in the time domain. (We could 
work in the frequency domain, using the stationary phase ap¬ 
proximation, but prefer to work in the time domain, since 
we consider the time domain versions of most approximants.) 
Nevertheless, if we take At = 0, so we just consider Atf, and 
also take the amplitudes of the two waveforms to be equal, 
then it is in fact straightforward to give a lower bound on the 
white noise match in terms of ||A(,c)|| 2 , if we further assume 
that we have a long enough observing time that we are able to 
discard rapidly varying contributions to the integrals. 

We thus hrst note that the white noise match between two 
waveforms hi and /i 2 is given by 


{hi,h2)2 - 


1 

2 

1 

2 



A^(v(t)) cos Aip(v(t))dt 
A^(v) cos A(p{u)t'(v)dv. 


(A5) 


Now, cos.r > 1 - x^/2 (which can be obtained by antidiffer¬ 
entiating cos X < 1 twice), so we have 


Mui{hi,h2) > 1 - 


\\Ah'\U\Aip\\2 
2||A o v\\l 


(A6) 


where we have used the inequality < II/IIOOII 0 II 2 , which 

can be obtained directly or from the Holder inequality. Here 
IIA^f'lloo and ||A o v \\2 are both hnite, since we have ymin > 0, 
so in this simple case, we can see that bounding ||Ai /:||2 to be 
close to zero implies that the matches should be close to unity. 
Note that we will consider norms over the full interval [0, ymaxl 
in the subsequent discussion, but the norms over the smaller 
interval we consider here are, of course, bounded from above 
by those over the larger interval. 

Thus, we note that we have ||Ay 2||2 < - 

0.64||Ai,c||oo. We can also obtain a bound on ||Ai,c )||2 in terms 
of ||Af|| 2 . Here we start from the same integrated-by-parts ex¬ 
pression for A(f we used above to obtain Eq. (A2) to now ob¬ 
tain 


M%A^\\, 


■ri' 

I 


[Af(y)]^ - 6 y^Af(y) I v^At(v)dv 
0 

2 x 


•f 


V At(v)dv 


dv 


^ cjmil 

X dv 


j|"“[6y^|Af(y)|||5||2||Af||2 + 9||5||2||Af||2] 


(yLx + 6||c||2lkll2 + 9y™axlkll2)l|Af|li 


(A7) 


Here we have dehned c(y) := y^ and s(v) (so ||c ||2 = 

fmax/ V7 and ||s ||2 = ymax/ VS) again used the triangle 
inequality (and triangle inequality for integrals) to obtain the 
first inequality, along with the Cauchy-Schwarz version of the 
Holder inequality, ||/g'||i < II/II2II0II2, from which we obtain 

r k(y)Af(y)|ufy < f |i(y)Af(y)|ufy = ||iAf||i < IkIbllAfIb 
Jo Jo 

(A8) 

(we also use the version with s —> c). We thus have 


MIIA^Ib < 




||Af||2-0.133||Af||2. 


(A9) 


Mw{hi,h 2 ) max 

tcVc 


(hi, 112)2 ^ { h \, h 2)2 
WhMMh ~ I|hi||2l|h2ll2’ 


(A4) 


In these I? inequalities, replacing the lower bound of y = 0 
by y = yg > 0 (as we actually do in the hts) would at most 
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just reduce the size of the constants in the hnal inequalities; 
The only change to the derivation is the different integration 
interval, since the basic integrated by parts expression does 
not change in this case, as discussed above. 


Appendix B: Prospects for improving more than one 
approximant with the same pPN coefflcients 

As a first example, we consider a function of the form 
Hiv) !F(f)/£'(y) for v e [0, ymax], where we know 'F and 
£ to some approximation. We want to know whether improv¬ 
ing the accuracy of H by adding additional terms to F (u) and 
&{v) will necessarily also improve the accuracy of other ap- 
proximants to consisting of 7'„['R] and 1/7’„[1/"??], where 
7’„[ ] denotes the nth order Maclaurin expansion (i.e., Tay¬ 
lor expansion about 0) of its argument. We thus note that 
if F(v) Aw“[l -I- f(c)] and &'{v) Bi^[\ + e(y)] with 

A, B e R and a, h 6 N constants, e and f real analytic, and 
|e(!;)|, |f(t;)| < 1 for y e [0, Umax], then one can obtain an ar¬ 
bitrarily good approximation to H from any of these approx- 
imants by using sufficiently high-order Maclaurin expansions 
of F and £ (and sufficiently large ri) since the Maclaurin se¬ 
ries for both F and £ converge to the functions and the series 
coming from expanding 1 /£' or IjF also converge, due to the 
assumption about the absolute values of e and f. 

Of course, this is not the situation we are actually in. For 
instance, we are actually hrst fitting to 7'„[1/"??], not We 
also know that we have ln"(u) terms in the expansion of both 
the energy and flux and that they both diverge at the light ring 
in the extreme mass-ratio limit (where v - 3“*^^ ^ 0.577), 
as discussed in, e.g., [17]. However, since we are considering 
t’max = t’isco, this divergence does not necessarily concern us, 
and we in fact find that the true energy flux (as obtained from 
NR, EOBNR, or black hole perturbation theory) is well within 
a factor of 2 of the Newtonian energy and flux over this range 
(see, e.g.. Fig. 4 in [42] in the equal-mass case and Fig. 1 in 
[43] in the extreme mass-ratio limit), so that we have |f(u)| < 1 
there. The binding energy in the test mass limit is also well- 
behaved up to uisco, we also have |e(u)| < 1 for v < uisco 
in that limit. (There does not appear to be any numerical rel¬ 
ativity data for the binding energy versus GW frequency for 
comparable-mass binaries against which to compare in the lit¬ 
erature.) 

Perhaps more importantly, we are not considering making 
Maclaurin expansions to arbitrary order; We are rather adding 
on a small number of additional terms to a known series (about 
V - 0, but not exactly a Maclaurin series due to the log¬ 
arithms) and taking n in the other approximants 7’„[‘R] and 
1/7’„[1/'R] to be consistent with the added terms. We are ob¬ 
taining the hrst six additional terms by (basically) minimizing 
the l} norm of the difference between the true value of 1 /T? 
and our 7’„[l/7?]. We then hx the hnal two additional terms 
by minimizing the difference between the true value of 7? and 
our (unexpanded) 7?. 

Actually, we minimize the antiderivatives of 1/7? and 
r„[l/‘^] for the hrst hts (the ones involving TaylorT2), but 
we ignore this slight complication here, since we are only in¬ 
terested in the accuracy of the antiderivatives themselves, as 
they are what gives the waveform, and controlling the function 
itself with an W norm is sufficient to control the antideriva¬ 
tive, since we are working on a compact set [cf Eq. (A8)]. 
Thus, our assumption about the smallness of the norm of 
1/7? - ^^[l/T?] is sufficient to guarantee that the norm of the 


antiderivatives will also be small. However, we cannot claim 
a priori that the ht of the antiderivatives will give us good 
agreement for the function itself (though this will likely be 
the case). 

As an example of the situations under which we can ex¬ 
pect to get an improvement of TaylorTl by fitting to im¬ 
prove TaylorT2, consider the case where we are just con¬ 
sidering ppN terms in the energy flux (which is the case 
at 4 and 4.5pPN, and can be taken to be the case at 5pPN 
when we are just looking at TaylorT2), so we are first fitting 
7’iv[£k„own/^£«pPN] to £eob/^eob and want to know if this 
will improve the fit of 7^<npPN/£known to Teob/Seob- We thus 
write F<„ppn = Tj^nown “ '?«pPN (where we introduce the minus 
sign to slightly simplify things later) and take N = 2n7 We 
thus want to bound ll^<«pPN/£known “ ^eob/SeobHp i" terms 
of 11 T2n [£]jnown /'^<npPN ] ~ ^eOB / ^OB lip, 117’2« [£i;nown ! ^nown ] ~ 
^known/^knownllp, and \\Fnp?ti\\p, in addition to known con¬ 
stants, where hats denote scaling by the Newtonian values (af¬ 
ter taking the derivative, for £') and we consider an arbitrary 
LP norm with p e [1, oo]. The Newtonian scaling is necessary 
for the first two norms to be finite, since &' jF goes as as 
u\0. 

We hrst note that we can write 


Tin 


^known 

F<np¥H 


- Tin 


a 


known 


Fkn 


+ FnpPN, 


(Bl) 


where Fn and £iv are the Newtonian Hux and binding energy, 
respectively. We thus have (adding a convenient zero) 


Tin 


^known 

TT/ipPN 


F’ 

^/eob 

Fpop 


- Tin 


^known 


Fkn. 




known 


Fkn 


■ FnpPli ■ 


^known 

Fk nown 


(B2) 

6/ 

^EOB 

T^eob 


Now, we can write the last three terms as (adding another con¬ 
venient zero) 


6 / 6 / 
^known^EOB 

Fk nown Feob 


Tkn own Feqb 


S' 


- 1 


EOB 


^^pPN ^ FeOB F^<npPN 




known 


F' F' 

*^EOB *^known 


(B3) 

Thus, applying the reverse triangle inequality twice, and not¬ 
ing that \\f/g\\p > ||/||p/||0||oc (which is just a rewriting of a 
case of the Holder inequality used earlier for p < oo and triv¬ 
ially true for p - oo), we hnd that 


Tin 


F' 

‘^EOB 




-^<«pPN. 

Teob 


■/TOB 

F<nppn 


^^\nown'^^OB ^ 

T^pPN 


F' 

‘^EOB 

^known 

P 

^EOB 

^known 

P- 


X 


Fk nown Feob 


F' F' 

*^known*^EOB 


Tin 

^known 

^known 




, -^known. 

-^known 


(B4) 


^ Here we use the subscript < «pPN to denote that we are including all the 
terms up to a given pPN order, which we just denoted by npPN previously, 
since we now reserve that subscript for just the «pPN terms. 
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Therefore, 




^known 

pf 

‘^EOB 


Tin 

^known 

1 


, ^ <«pPN . 


F’ 

‘-^EOB 


+ 

P 


^known 

^known 

^ In 

, ^^^laiown. 




-H 

^^nown^^EOB 

- 1 

^^^pPN 

^known^EOB 

OO 

^EOB 


^known 


(B5) 


This inequality seemingly implies that in order for us 
to have our fit for TaylorT2 also improve TaylorTl (or 
TaylorT4), we should also have the higher PN expan¬ 
sion of fi^nown/'^known ^gfc^ing Well with the unexpanded 
^Lown/^kiiown and also a small (and we indeed saw that 
smaller pPN coefficients helped the TaylorT2 fits to also im¬ 


prove the agreement of TaylorTl and—to a lesser extent— 
TaylorT4). However, this is just an upper bound, and is likely 
not sharp, particularly since we have used the weaker version 
of the reverse triangle inequality without the absolute value 
(i.e., the first and last terms in ||x- i/|| > |||x|| - ||^||| > ||jc|| - ||i/|| 
instead of the first and second terms). Indeed, at 4pPN, the 
right-hand side of Eq. (B5) is ~ 50-100 times larger than the 
left-hand side if we take p - 1, 2, or cm. The difference be¬ 
tween the two sides is not so pronounced (< 50) if we con¬ 
sider the case where we include the Newtonian scaling on all 
terms, where the inequality also holds. Thus, one cannot use 
this equality to say that a small !F),pPN is necessary for the Tay- 
lorT2 fit to also improve TaylorTl, though this is certainly 
a sufficient condition (along with having appropriately small 
values for the other two norms). Indeed, in the 4pPN exam¬ 
ples mentioned above, the contribution from the final term is 
at most a quarter of the total, and often much less. 
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